Libraries

library("corrplot")
## corrplot 0.95 loaded
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
dsMERGIATO <- read.csv("/Users/leomarcellopoli/Documents/Nonparametric/Project/Non-parametric-Statistics-2025-2026/dataset/dataset_MERGIATO", header = TRUE, sep = ",")

Import dataset

set.seed(123) # Per riproducibilità
indici_random <- sample(1:nrow(dsMERGIATO), 10000)
dsridotto <- dsMERGIATO[indici_random, -1]      # e senza date
dsclorofilla <- dsridotto$Chl
dscovariate <- dsridotto[, -which(colnames(dsridotto) == "Chl")]

Correlation plot

corrplot(cor(dsridotto), method = "circle")

Scatterplot

for (i in 1:ncol(dscovariate)) {
  
  # A. Estrazione dinamica del nome della colonna X
  nome_covariata <- colnames(dscovariate)[i]
  
  # B. Plot con xlab dinamico
  # plot(X, Y, ...)
  plot(dscovariate[,i], 
       dsclorofilla, 
       pch = 16, 
       cex = 0.5, # Punti piccoli
       col = rgb(0.7, 0, 0, 0.4), # Colore rosso trasparente
       ylab = "Clorofilla (Chl)",
       xlab = nome_covariata, # <--- QUI STA LA CORREZIONE
       main = paste("Chl vs", nome_covariata)
  ) 
}

PCA su tutte le covariate

# vars_to_keep <- c("Temp", "Salinity", "MLD", "uo", "vo", "Current_Speed", "Solar_Flux", "Heat_Flux")
vars_to_keep <- c("Temp", "Salinity", "MLD", "Current_Speed", "Solar_Flux", "Heat_Flux")
data_for_pca <- dsridotto %>%
  select(all_of(vars_to_keep))

# Boxplot dati grezzi (Vedrai che MLD o Flux schiacciano tutto il resto)
par(mfrow = c(1, 2))
boxplot(data_for_pca, las = 2, col = 'gold', main = "Dati Grezzi (Scale diverse!)")

# Boxplot dati standardizzati (Ora sono confrontabili)
data_scaled <- scale(data_for_pca)
boxplot(data_scaled, las = 2, col = 'orange', main = "Dati Standardizzati")

par(mfrow = c(1, 1)) # Ripristina layout

pca_result <- prcomp(data_scaled, center = TRUE, scale. = TRUE)
pca_result
## Standard deviations (1, .., p=6):
## [1] 1.6666368 1.1043476 0.8780764 0.8080498 0.6368531 0.4161652
## 
## Rotation (n x k) = (6 x 6):
##                      PC1        PC2         PC3        PC4         PC5
## Temp           0.3760473 -0.3364457 -0.57528592 -0.5439980  0.03135252
## Salinity      -0.3360247 -0.4925235 -0.42775929  0.5016859 -0.45805260
## MLD           -0.4795047 -0.2364082 -0.12887436  0.0839011  0.81980387
## Current_Speed  0.1792844  0.6577851 -0.59243644  0.3962656  0.14928184
## Solar_Flux     0.5066973 -0.3296838 -0.08306083  0.2118105  0.28400341
## Heat_Flux      0.4763045 -0.2167104  0.33403512  0.4934045  0.11909427
##                        PC6
## Temp          -0.342814895
## Salinity      -0.005388042
## MLD           -0.135867614
## Current_Speed -0.069881806
## Solar_Flux     0.708622402
## Heat_Flux     -0.597461073
layout(matrix(c(1, 2), 1, 2))
# Plot degli autovalori (Varianza di ogni PC)
{
plot(pca_result, type = "lines", main = "Scree Plot")
abline(h = 1, col = "red", lty = 2) # Criterio di Kaiser: tieni PC con autovalore > 1
}

cum_var <- cumsum(pca_result$sdev^2) / sum(pca_result$sdev^2)
{
plot(cum_var, type = "b", xlab = "Numero di Componenti", 
     ylab = "Varianza Cumulata", ylim = c(0, 1), main = "Varianza Spiegata Totale")
abline(h = 0.90, col = "blue", lty = 2) # Linea del 90%
}

layout(1) # Reset

# --- 5. Interpretazione: I Loadings (Il "DNA" delle componenti) ---

# Carichiamo i coefficienti dalla variabile 'rotation'
# NOTA: Usiamo 'rotation' invece di 'loadings' per prcomp()
loadings_mat <- pca_result$rotation 
# Visualizziamo i pesi delle prime 3 componenti
par(mfrow = c(3, 1), mar = c(2, 4, 2, 1))
# Ciclismo sulle prime 3 componenti principali (PC1, PC2, PC3)
for(i in 1:min(3, ncol(loadings_mat))) {
  # Barplot: Passiamo la colonna 'i' della matrice di rotazione
  barplot(loadings_mat[, i], 
          ylim = c(-1, 1), 
          main = paste("Loadings PC", i), 
          col = "lightblue",
          las = 2 # Ruota le etichette dell'asse X (variabili) per renderle leggibili
  )
  # Aggiunge la linea orizzontale a zero per facilitare la lettura dei segni (+/-)
  abline(h = 0)
}

# Ripristina il layout grafico
par(mfrow = c(1, 1))

# Mostra sia le osservazioni (punti) che le variabili (frecce)
# - Frecce vicine = Variabili molto correlate
# - Frecce lunghe = Variabili importanti per queste componenti
biplot(pca_result, cex = 0.7, col = c("gray", "red"), 
       main = "Biplot: Relazioni tra Variabili")

# --- 7. (Opzionale) Salva gli Scores nel dataset originale ---
# Aggiungi le nuove "Super-Variabili" al tuo dataset per usarle nella regressione
scores <- as.data.frame(pca_result$scores)